
make_scale_coefs <- function(.df, .dep_var, ...){
  
  #dep vars scaling
  .dep_var = enquo(.dep_var)
  .df =
    .df %>%
    mutate(outcome = !!.dep_var) %>%
    mutate(scaled_dv_outcome = scale(outcome),
           scaled_iv_pop = scale(log10(pop)))
  
  #independent vars scaling
  .ind_vars = enquos(...)
  for(i in .ind_vars){
    .df =
      .df %>%
      mutate(ind_var = scale(!!i))
    
    names(.df)[names(.df) == "ind_var"] = paste0("scaled_iv_", quo_name(i))
    
  }
  
  
  form_reg <-
    as.formula(
      paste0("scaled_dv_outcome ~ scaled_iv_pop +",
             tidyselect::vars_select(names(.df), starts_with('scaled_iv_', ignore.case = TRUE)) %>%
               paste0(., collapse = " + ")))
  
  lm(formula = form_reg, data = .df) %>%
    tidy() %>%
    return()
  
}


doc_mod_indo <- make_scale_coefs(.df = indonesia_data, .dep_var = prop_doctors, ... = pop_growth, density, div_index, avg_income, elec_comp_index) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Indonesia") %>% mutate(outcome = "docto")
teacher_mod_indo <- make_scale_coefs(.df = indonesia_data, .dep_var = teachers_prop, ... = pop_growth, density, div_index, avg_income, elec_comp_index) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Indonesia") %>% mutate(outcome = "teacher")
comm_mod_indo <- make_scale_coefs(.df = indonesia_data, .dep_var = prop_comm_health, ... = pop_growth, density, div_index, avg_income, elec_comp_index) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Indonesia") %>% mutate(outcome = "health")
schools_mod_indo <- make_scale_coefs(.df = indonesia_data, .dep_var = schools_prop, ... = pop_growth, density, div_index, avg_income, elec_comp_index) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Indonesia") %>% mutate(outcome = "schools")
elec_mod_indo <- make_scale_coefs(.df = indonesia_data, .dep_var = prop_govt_electricity, ... = pop_growth, density, div_index, avg_income, elec_comp_index) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Indonesia") %>% mutate(outcome = "elect")
water_mod_indo <- make_scale_coefs(.df = indonesia_data, .dep_var = prop_water, ... = pop_growth, density, div_index, avg_income, elec_comp_index) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Indonesia") %>% mutate(outcome = "water")

#subset on urban municipalities
brazil_data_models <- brazil_data[brazil_data$urb == 1,]

doc_mod_brazil <- make_scale_coefs(.df = brazil_data_models, .dep_var = prop_doctos, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "docto") 
teacher_mod_brazil <- make_scale_coefs(.df = brazil_data_models, .dep_var = teachers_prop, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "teacher") 
comm_mod_brazil <- make_scale_coefs(.df = brazil_data_models, .dep_var = prop_comm_health, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "health") 
schools_mod_brazil <- make_scale_coefs(.df = brazil_data_models, .dep_var = schools_prop, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "schools") 
elec_mod_brazil <- make_scale_coefs(.df = brazil_data_models, .dep_var = electr_share_2010, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "elect") 
water_mod_brazil <- make_scale_coefs(.df = brazil_data_models, .dep_var = piped_water_share_2010, ... = pop_growth, density, div_index, avg_salary_mo_2010, ill_rate_2010, elec_index, align) %>%
  filter(term == "scaled_iv_pop") %>% mutate(country = "Brazil") %>% mutate(outcome = "water") 


coefficient_plot <-
  bind_rows(comm_mod_indo, schools_mod_indo, doc_mod_indo, teacher_mod_indo, elec_mod_indo, water_mod_indo,
            comm_mod_brazil, schools_mod_brazil,doc_mod_brazil, teacher_mod_brazil, elec_mod_brazil, water_mod_brazil) %>%
  mutate(label = case_when(str_detect(outcome, "schools") ~ "Schools (per 1000)",
                           str_detect(outcome, "teacher") ~ "Teachers (per 1000)",
                           str_detect(outcome, "water") ~ "Piped Water (%)",
                           str_detect(outcome, "elect") ~ "Electricity Access (%)",
                           str_detect(outcome, "docto") ~ "Doctors (per 1000)",
                           str_detect(outcome, "health") ~ "Health Centers (per 1000)",
                           TRUE ~ NA_character_)) %>%
  mutate(rank = case_when(str_detect(outcome, "schools") ~ "2",
                          str_detect(outcome, "water") ~ "9",
                          str_detect(outcome, "elect") ~ "8",
                          str_detect(outcome, "teacher") ~ "6",
                          str_detect(outcome, "docto") ~ "5",
                          str_detect(outcome, "health") ~ "3",
                          TRUE ~ NA_character_)) %>%
  dplyr::select(estimate, std.error, label, country, rank) %>%
  add_row(estimate = NA_real_, std.error = NA_real_, label = "Divisible Public Goods:                             ", country = "Brazil", rank = "1") %>%
  add_row(estimate = NA_real_, std.error = NA_real_, label = "Personnel:                                          ", country = "Brazil", rank = "4") %>%
  add_row(estimate = NA_real_, std.error = NA_real_, label = "Indivisible Public Goods:                           ", country = "Brazil", rank = "7") %>%
  mutate(rank = as.numeric(rank)) %>%
  ggplot(aes(x=estimate, y = reorder(label, -rank))) +
  geom_point() +
  geom_errorbarh(aes(xmin=estimate-1.96*std.error, xmax = estimate+1.96*std.error), height = 0.1) +
  geom_vline(xintercept = 0, linetype = "dotted") +
  facet_wrap(country~.) +
  theme_bw() +
  scale_colour_grey() +
  theme(panel.grid.minor = element_blank(), 
        panel.grid.major.x = element_blank(),
        axis.line.y.left = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "bottom",
        strip.background = element_blank(),
        legend.title = element_blank(),
        axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        axis.title.y = element_blank()) +
  xlab("Standardized Coefficient")


ggsave("./_4_outputs/figure_3.pdf", plot = coefficient_plot, width = 8, height = 4)


